clear all
set more off
set mem 10000000
set matsize 10000
version 15

*****************************************************************************************
*** RDROBUST first stage, nighttime brightness: splits on hours of rural power supply ***
*****************************************************************************************

** Set file paths
do "$path_code/paths.do"

** Set graph scheme
cd "$path/code/analyze"
set scheme fb, perm

****************************************************************** 
****************************************************************** 

{
use "$panel/panel_dataset_full.dta", clear
gen in_fs_sample = vplan4<11 & corr_state==1 & sample==1 & sample_h==1 & sing_h==1 & pop_non_zero==1 

	// Summarize hours of power to electrified villages, by district
egen temp1 = mean(vdp_pwr_h_all_avg) if vdp_pwr_h_all_avg_11>0, by(st_code dt_code)
egen HRS_all_wide = mode(temp1), by(st_code dt_code)
egen temp2 = mean(vdp_pwr_h_dom_avg) if vdp_pwr_h_dom_avg_11>0, by(st_code dt_code)
egen HRS_dom_wide = mode(temp2), by(st_code dt_code)
egen temp3 = mean(vdp_pwr_h_all_avg) if vdp_pwr_h_all_avg_11>0 & inrange(tot_p,150,450), by(st_code dt_code)
egen HRS_all_rd = mode(temp3), by(st_code dt_code)
egen temp4 = mean(vdp_pwr_h_dom_avg) if vdp_pwr_h_dom_avg_11>0  & inrange(tot_p,150,450), by(st_code dt_code)
egen HRS_dom_rd = mode(temp4), by(st_code dt_code)

twoway scatter HRS_all_wide HRS_dom_wide if in_fs_sample==1
twoway scatter HRS_all_rd HRS_dom_rd if in_fs_sample==1
twoway scatter HRS_all_wide HRS_all_rd if in_fs_sample==1
twoway scatter HRS_dom_wide HRS_dom_rd if in_fs_sample==1

	// High-power-quality states
gen hpq_st = 0  // based on state-level power deficits as reported in the Load Generation Balance Report: 2011-12 http://www.cea.nic.in/reports/annual/lgbr/lgbr-2011.pdf
replace hpq_st = 1 if state=="SIKKIM" // 90.5
replace hpq_st = 1 if state=="CHHATTISGARH" // 17.3
replace hpq_st = 1 if state=="TRIPURA" // 13.0
replace hpq_st = 1 if state=="HIMACHAL PRADESH" // 7.1
replace hpq_st = 1 if state=="KARNATAKA" // 4.8
replace hpq_st = 1 if state=="MIZORAM" // 4.2
replace hpq_st = 1 if state=="WEST BENGAL" // 0.0
replace hpq_st = 1 if state=="ASSAM" // -0.8
replace hpq_st = 1 if state=="MANIPUR" // -0.9
replace hpq_st = 1 if state=="ARUNACHAL PRADESH" // -1.1
replace hpq_st = 1 if state=="GUJARAT" // -1.6
replace hpq_st = 1 if state=="MEGHALAYA" // -2.7
replace hpq_st = 1 if state=="HARYANA" // -6.0
replace hpq_st = 1 if state=="RAJASTHAN" // -7.0
tab state hpq_st
tab state hpq_st if inrange(tot_p,100,500)
tab state hpq_st if in_fs_sample==1
tab state hpq_st if inrange(tot_p,100,500) & in_fs_sample==1

tabstat HRS_all_wide if inrange(tot_p,150,450) & hpq_st==1 & in_fs_sample==1, by(state) s(p10 p25 p50 p75 p90)
tabstat HRS_all_wide if inrange(tot_p,150,450) & hpq_st==0 & in_fs_sample==1, by(state) s(p10 p25 p50 p75 p90)

count if inrange(tot_p,150,450) & hpq_st==1 & in_fs_sample==1 & HRS_all_wide>10 & HRS_all_wide!=.
count if inrange(tot_p,150,450) & hpq_st==1 & in_fs_sample==1 & HRS_all_wide<=10 & HRS_all_wide!=.

count if inrange(tot_p,150,450) & hpq_st==0 & in_fs_sample==1 & HRS_all_wide>10 & HRS_all_wide!=.
count if inrange(tot_p,150,450) & hpq_st==0 & in_fs_sample==1 & HRS_all_wide<=10 & HRS_all_wide!=.

count if inrange(tot_p,150,450) & in_fs_sample==1 & HRS_all_wide>=10 & HRS_all_wide!=.
count if inrange(tot_p,150,450) & in_fs_sample==1 & HRS_all_wide<10 & HRS_all_wide!=.

	// Keep villages in RD sample
keep if in_fs_sample==1
	
	// Create state FEs (since RD robust doesn't let you pass them through)
drop if st_code==32 // Kerala, only 3 villages
tab st_code, gen(STFE)	
drop STFE1 // to avoid collinearity

	// Create district FEs (since RD robust doesn't let you pass them through)
gen stdtFE = stdt
tab stdtFE state     if inlist(stdtFE,62,63,64,70,91,180,181,202,308,309,320,490,493,494,495,504,505,508)	
replace stdtFE = 999 if inlist(stdtFE,62,63,64,70,91,180,181,202,308,309,320,490,493,494,495,504,505,508)	
	// one catch-all district FE for districts with so few in-sample villages that they break rdrobust
tab stdtFE, gen(DTFE)	
drop DTFE1 // to avoid collinearity

	// Create block groups (for clustering)
egen stdtbk = group(stdt bk_code)

	// Create lights-difference variable, to identify crazy outliers
gen lights_diff = abs(lights_max2011_hat - lights_max2001_hat)


  // Create variables to store regresson results
gen fs_step = .
gen hrs_step = .
gen yvar = ""
gen ifs = ""
gen hrs_if = ""
gen control = ""
gen fe = ""
gen kernel = ""
gen bwmethod = ""
gen vce = ""
gen polynomial_order = .
gen beta_conv = .
gen beta_robust = .
gen se_conv = .
gen se_robust = .
gen pval_conv = .
gen pval_robust = .
gen lci_conv = .
gen uci_conv = .
gen lci_robust = .
gen uci_robust = .
gen bw_lo = .
gen bw_hi = .
gen nobs_orig = .
gen nobs_left = .
gen nobs_right = .
gen ndist = . 
gen ymean = .
gen ftag = ""
gen hrstag = ""
local folder = "RDROBUST plots fs lights HRS"
local h_wide = 200


	// 1.1 Lights brightness regressions, preferred
local fs_step = 1
local ifs = "pop_mismatch20==0 & lights_diff<20"
local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
local fe = "STFE*"
local kernel = "tri"
local bwmethod = "mserd"
local vce = ""
local poly = 1
local row = 0

foreach hrs in 1 2 {
	
	if `hrs'==1 {
		local hrs_if = " & HRS_all_wide>=10 & HRS_all_wide!=."
		local hrstag = "over10hrs"
	}
	else if `hrs'==2 {
		local hrs_if = " & HRS_all_wide<10 & HRS_all_wide!=."
		local hrstag = "under10hrs"
    }
   
	foreach y in 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 {

		// Define outcome variable
		local yvar = "lights_max`y'_hat"
		local title = "Preferred specification: `y' Brightness"
		local ftag = "preferred"

		// Remove controls used to project outcome variable
		local controls = "`CONTROLS'"
		if regexm("`yvar'","2007")==1 {
			local controls = subinstr("`controls'"," lights_max2005","",1)
		}
		else if regexm("`yvar'","2006")==1 {
			local controls = subinstr("`controls'"," lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2005")==1 {
			local controls = subinstr("`controls'"," lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2004")==1 {
			local controls = subinstr("`controls'"," lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2003")==1 {
			local controls = subinstr("`controls'"," lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2002")==1 {
			local controls = subinstr("`controls'"," lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		
		// Run first-stage regresssion
		rdrobust `yvar' tot_p if `ifs' `hrs_if', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

		// Generate in-sample indicator and store bandwidths
		cap drop temp_in_reg
		qui gen temp_in_reg = `ifs' `hrs_if' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r))
		qui count if temp_in_reg
		assert r(N)==(e(N_h_l)+e(N_h_r))
		local h_l = e(h_l)
		local h_r = e(h_r)

		// Store results
		local row = `row' + 1
		qui replace fs_step = `fs_step' in `row'
		qui replace hrs_step = `hrs' in `row'
		qui replace yvar = "`yvar'" in `row'
		qui replace ifs = "`ifs'" in `row'
		qui replace hrs_if = "`hrs_if'" in `row'
		qui replace control = "`controls'" in `row'
		qui replace fe = "`fe'" in `row'
		qui replace kernel = e(kernel) in `row'
		qui replace bwmethod = e(bwselect) in `row'
		qui replace vce = "`vce'" in `row'
		qui replace polynomial_order = e(p) in `row'
		qui replace beta_conv = e(tau_cl) in `row'
		qui replace beta_robust = e(tau_bc) in `row'
		qui replace se_conv = e(se_tau_cl) in `row'
		qui replace se_robust = e(se_tau_rb) in `row'
		qui replace pval_conv = e(pv_cl) in `row'
		qui replace pval_robust = e(pv_rb) in `row'
		qui replace lci_conv = e(ci_l_cl) in `row'
		qui replace uci_conv = e(ci_r_cl) in `row'
		qui replace lci_robust = e(ci_l_rb) in `row'
		qui replace uci_robust = e(ci_r_rb) in `row'
		qui replace bw_lo = e(h_l) in `row'
		qui replace bw_hi = e(h_r) in `row'
		qui replace nobs_orig = e(N) in `row'
		qui replace nobs_left = e(N_h_l) in `row'
		qui replace nobs_right = e(N_h_r) in `row'
		qui unique stdt if temp_in_reg==1
		qui replace ndist = r(unique) in `row'
		qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
		qui replace ymean = r(mean) in `row'
		qui replace ftag = "`ftag'" in `row'
		qui replace hrstag = "`hrstag'" in `row'

		// Residualize `yvar' for regression sample
		qui reg `yvar' `controls' `fe' if temp_in_reg==1
		cap drop y_resid_sample
		predict y_resid_sample, residuals
		
		// Residualize `yvar' for [100,500] bandwidth
		qui reg `yvar' `controls' `fe' if `ifs' `hrs_if' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
		cap drop y_resid_full
		predict y_resid_full, residuals
		
		// RD plots for regression sample (variable bandwidth) 
		foreach bins in 10 20 40 {
			if `h_l'<=75 {
				local xlab = "250 275 300 325 350"
			}
			else if inrange(`h_l',76,125) {
				local xlab = "200 250 300 350 400"
			}
			else if inrange(`h_l',126,175) {
				local xlab = "150 200 250 300 350 400 450"
			}
			else {
				local xlab = ""
			}
			rdplot y_resid_sample tot_p if temp_in_reg==1, ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`y' brightness residuals", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) ///
				xlabel(`xlab', labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/lights_fs_`ftag'_`hrstag'_`y'_inreg_`bins'.pdf", replace
		}

		// RD plots for regression sample (constant bandwidth) 
		foreach bins in 10 20 40 {
			if `bins'==40 {
				local ymarks = "-.4 -.2 0 .2 .4"
			}
			else {
				local ymarks = "-.3 -.2 -.1 0 .1 .2 .3"
			}
			rdplot y_resid_full tot_p if `ifs' `hrs_if' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`y' brightness residuals", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(`ymarks',nogrid angle(0) labsize(medlarge)) ///
				xlabel(, labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/lights_fs_`ftag'_`hrstag'_`y'_wide_`bins'.pdf", replace
			local ymarks = ""
		}			
	}
}

	// 1.7 Sensitivity to different bandwidth calculations
local fs_step = 7
local ifs = "pop_mismatch20==0 & lights_diff<20"
local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
local fe = "STFE*"
local kernel = "tri"
local vce = ""
local poly = 1

foreach hrs in 1 2 {
	
	if `hrs'==1 {
		local hrs_if = " & HRS_all_wide>=10 & HRS_all_wide!=."
		local hrstag = "over10hrs"
	}
	else if `hrs'==2 {
		local hrs_if = " & HRS_all_wide<10 & HRS_all_wide!=."
		local hrstag = "under10hrs"
    }
   
	foreach y in 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 {

		foreach reg in 1 2 3 4 5 {

			// Define outcome variable
			local yvar = "lights_max`y'_hat"

			// Define bandwidth method
			if `reg'==1 {
				local bwmethod = "msetwo"
				local title = "`y' Brightness, Bandwidth Method `bwmethod'"
				local ftag = "BW`bwmethod'"
			}
			if `reg'==2 {
				local bwmethod = "msesum"
				local title = "`y' Brightness, Bandwidth Method `bwmethod'"
				local ftag = "BW`bwmethod'"
			}
			if `reg'==3 {
				local bwmethod = "cerrd"
				local title = "`y' Brightness, Bandwidth Method `bwmethod'"
				local ftag = "BW`bwmethod'"
			}
			if `reg'==4 {
				local bwmethod = "certwo"
				local title = "`y' Brightness, Bandwidth Method `bwmethod'"
				local ftag = "BW`bwmethod'"
			}
			if `reg'==5 {
				local bwmethod = "cersum"
				local title = "`y' Brightness, Bandwidth Method `bwmethod'"
				local ftag = "BW`bwmethod'"
			}

			
			// Remove controls used to project outcome variable
			local controls = "`CONTROLS'"
			if regexm("`yvar'","2007")==1 {
				local controls = subinstr("`controls'"," lights_max2005","",1)
			}
			else if regexm("`yvar'","2006")==1 {
				local controls = subinstr("`controls'"," lights_max2004 lights_max2005","",1)
			}
			else if regexm("`yvar'","2005")==1 {
				local controls = subinstr("`controls'"," lights_max2003 lights_max2004 lights_max2005","",1)
			}
			else if regexm("`yvar'","2004")==1 {
				local controls = subinstr("`controls'"," lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
			}
			else if regexm("`yvar'","2003")==1 {
				local controls = subinstr("`controls'"," lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
			}
			else if regexm("`yvar'","2002")==1 {
				local controls = subinstr("`controls'"," lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
			}
			
			// Run first-stage regresssion
			rdrobust `yvar' tot_p if `ifs' `hrs_if', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

			// Generate in-sample indicator and store bandwidths
			cap drop temp_in_reg
			qui gen temp_in_reg = `ifs' `hrs_if' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r))
			qui count if temp_in_reg
			assert r(N)==(e(N_h_l)+e(N_h_r))
			local h_l = e(h_l)
			local h_r = e(h_r)
			
			// Store results
			local row = `row' + 1
			qui replace fs_step = `fs_step' in `row'
			qui replace hrs_step = `hrs' in `row'
			qui replace yvar = "`yvar'" in `row'
			qui replace ifs = "`ifs'" in `row'
			qui replace hrs_if = "`hrs_if'" in `row'
			qui replace control = "`controls'" in `row'
			qui replace fe = "`fe'" in `row'
			qui replace kernel = e(kernel) in `row'
			qui replace bwmethod = e(bwselect) in `row'
			qui replace vce = "`vce'" in `row'
			qui replace polynomial_order = e(p) in `row'
			qui replace beta_conv = e(tau_cl) in `row'
			qui replace beta_robust = e(tau_bc) in `row'
			qui replace se_conv = e(se_tau_cl) in `row'
			qui replace se_robust = e(se_tau_rb) in `row'
			qui replace pval_conv = e(pv_cl) in `row'
			qui replace pval_robust = e(pv_rb) in `row'
			qui replace lci_conv = e(ci_l_cl) in `row'
			qui replace uci_conv = e(ci_r_cl) in `row'
			qui replace lci_robust = e(ci_l_rb) in `row'
			qui replace uci_robust = e(ci_r_rb) in `row'
			qui replace bw_lo = e(h_l) in `row'
			qui replace bw_hi = e(h_r) in `row'
			qui replace nobs_orig = e(N) in `row'
			qui replace nobs_left = e(N_h_l) in `row'
			qui replace nobs_right = e(N_h_r) in `row'
			qui unique stdt if temp_in_reg==1
			qui replace ndist = r(unique) in `row'
			qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
			qui replace ymean = r(mean) in `row'
			qui replace ftag = "`ftag'" in `row'
			qui replace hrstag = "`hrstag'" in `row'

			// Residualize `yvar' for regression sample
			qui reg `yvar' `controls' `fe' if temp_in_reg==1
			cap drop y_resid_sample
			predict y_resid_sample, residuals
			
			// RD plots for regression sample (variable bandwidth) 
			foreach bins in 10 20 40 {
				if `h_l'<=75 {
					local xlab = "250 275 300 325 350"
				}
				else if inrange(`h_l',76,125) {
					local xlab = "200 250 300 350 400"
				}
				else if inrange(`h_l',126,175) {
					local xlab = "150 200 250 300 350 400 450"
				}
				else {
					local xlab = ""
				}
				rdplot y_resid_sample tot_p if temp_in_reg==1, ///
					c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
					graph_options( ///
					title("`title'", color(black) size(large)) ///
					ytitle("`y' brightness residuals", size(medlarge)) ///
					xtitle("2001 village population", size(medlarge)) ///
					ylabel(,nogrid angle(0) labsize(medlarge)) ///
					xlabel(`xlab', labsize(medlarge)) ///
					graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
					legend(off))
				graph export "$results/`folder'/lights_fs_`ftag'_`hrstag'_`y'_inreg_`bins'.pdf", replace
			}
		}		
	}
}


	// 1.8 Sensitivity to 2nd-order polynomial and differet kernels
local fs_step = 8
local ifs = "pop_mismatch20==0 & lights_diff<20"
local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
local fe = "STFE*"
local bwmethod = "mserd"
local vce = ""

foreach hrs in 1 2 {
	
	if `hrs'==1 {
		local hrs_if = " & HRS_all_wide>=10 & HRS_all_wide!=."
		local hrstag = "over10hrs"
	}
	else if `hrs'==2 {
		local hrs_if = " & HRS_all_wide<10 & HRS_all_wide!=."
		local hrstag = "under10hrs"
    }
   
	foreach y in 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 {

		foreach reg in 2 4 {

			// Define outcome variable
			local yvar = "lights_max`y'_hat"

			// Define kernel and polynomial combo
			if `reg'==1 {
				local kernel = "tri"
				local poly = 2
				local title = "`y' Brightness, Triangular Kernel, Quadratic"
				local ftag = "K`kernel'P`poly'"
			}
			if `reg'==2 {
				local kernel = "epa"
				local poly = 1
				local title = "`y' Brightness, Epanechnikov Kernel"
				local ftag = "K`kernel'P`poly'"
			}
			if `reg'==3 {
				local kernel = "epa"
				local poly = 2
				local title = "`y' Brightness, Epanechnikov Kernel, Quadratic"
				local ftag = "K`kernel'P`poly'"
			}
			if `reg'==4 {
				local kernel = "uni"
				local poly = 1
				local title = "`y' Brightness, Uniform Kernel"
				local ftag = "K`kernel'P`poly'"
			}
			if `reg'==5 {
				local kernel = "uni"
				local poly = 2
				local title = "`y' Brightness, Uniform Kernel, Quadratic"
				local ftag = "K`kernel'P`poly'"
			}

			
			// Remove controls used to project outcome variable
			local controls = "`CONTROLS'"
			if regexm("`yvar'","2007")==1 {
				local controls = subinstr("`controls'"," lights_max2005","",1)
			}
			else if regexm("`yvar'","2006")==1 {
				local controls = subinstr("`controls'"," lights_max2004 lights_max2005","",1)
			}
			else if regexm("`yvar'","2005")==1 {
				local controls = subinstr("`controls'"," lights_max2003 lights_max2004 lights_max2005","",1)
			}
			else if regexm("`yvar'","2004")==1 {
				local controls = subinstr("`controls'"," lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
			}
			else if regexm("`yvar'","2003")==1 {
				local controls = subinstr("`controls'"," lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
			}
			else if regexm("`yvar'","2002")==1 {
				local controls = subinstr("`controls'"," lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
			}
			
			// Run first-stage regresssion
			rdrobust `yvar' tot_p if `ifs' `hrs_if', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

			// Generate in-sample indicator and store bandwidths
			cap drop temp_in_reg
			qui gen temp_in_reg = `ifs' `hrs_if' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r))
			qui count if temp_in_reg
			assert r(N)==(e(N_h_l)+e(N_h_r))
			local h_l = e(h_l)
			local h_r = e(h_r)
			
			// Store results
			local row = `row' + 1
			qui replace fs_step = `fs_step' in `row'
			qui replace hrs_step = `hrs' in `row'
			qui replace yvar = "`yvar'" in `row'
			qui replace ifs = "`ifs'" in `row'
			qui replace hrs_if = "`hrs_if'" in `row'
			qui replace control = "`controls'" in `row'
			qui replace fe = "`fe'" in `row'
			qui replace kernel = e(kernel) in `row'
			qui replace bwmethod = e(bwselect) in `row'
			qui replace vce = "`vce'" in `row'
			qui replace polynomial_order = e(p) in `row'
			qui replace beta_conv = e(tau_cl) in `row'
			qui replace beta_robust = e(tau_bc) in `row'
			qui replace se_conv = e(se_tau_cl) in `row'
			qui replace se_robust = e(se_tau_rb) in `row'
			qui replace pval_conv = e(pv_cl) in `row'
			qui replace pval_robust = e(pv_rb) in `row'
			qui replace lci_conv = e(ci_l_cl) in `row'
			qui replace uci_conv = e(ci_r_cl) in `row'
			qui replace lci_robust = e(ci_l_rb) in `row'
			qui replace uci_robust = e(ci_r_rb) in `row'
			qui replace bw_lo = e(h_l) in `row'
			qui replace bw_hi = e(h_r) in `row'
			qui replace nobs_orig = e(N) in `row'
			qui replace nobs_left = e(N_h_l) in `row'
			qui replace nobs_right = e(N_h_r) in `row'
			qui unique stdt if temp_in_reg==1
			qui replace ndist = r(unique) in `row'
			qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
			qui replace ymean = r(mean) in `row'
			qui replace ftag = "`ftag'" in `row'
			qui replace hrstag = "`hrstag'" in `row'

			// Residualize `yvar' for regression sample
			qui reg `yvar' `controls' `fe' if temp_in_reg==1
			cap drop y_resid_sample
			predict y_resid_sample, residuals
			
			// Residualize `yvar' for [100,500] bandwidth
			qui reg `yvar' `controls' `fe' if `ifs' `hrs_if' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
			cap drop y_resid_full
			predict y_resid_full, residuals
			
			// RD plots for regression sample (variable bandwidth) 
			foreach bins in 10 20 40 {
				if `h_l'<=75 {
					local xlab = "250 275 300 325 350"
				}
				else if inrange(`h_l',76,125) {
					local xlab = "200 250 300 350 400"
				}
				else if inrange(`h_l',126,175) {
					local xlab = "150 200 250 300 350 400 450"
				}
				else {
					local xlab = ""
				}
				rdplot y_resid_sample tot_p if temp_in_reg==1, ///
					c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
					graph_options( ///
					title("`title'", color(black) size(large)) ///
					ytitle("`y' brightness residuals", size(medlarge)) ///
					xtitle("2001 village population", size(medlarge)) ///
					ylabel(,nogrid angle(0) labsize(medlarge)) ///
					xlabel(`xlab', labsize(medlarge)) ///
					graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
					legend(off))
				graph export "$results/`folder'/lights_fs_`ftag'_`hrstag'_`y'_inreg_`bins'.pdf", replace
			}

			// RD plots for regression sample (constant bandwidth) 
			foreach bins in 10 20 40 {
				rdplot y_resid_full tot_p if `ifs' `hrs_if' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
					c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
					graph_options( ///
					title("`title'", color(black) size(large)) ///
					ytitle("`y' brightness residuals", size(medlarge)) ///
					xtitle("2001 village population", size(medlarge)) ///
					ylabel(,nogrid angle(0) labsize(medlarge)) ///
					xlabel(, labsize(medlarge)) ///
					graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
					legend(off))
				graph export "$results/`folder'/lights_fs_`ftag'_`hrstag'_`y'_wide_`bins'.pdf", replace
			}
		}		
	}
}


	// Save results
keep fs_step-hrstag
dropmiss, obs force
compress
save "$results/RDROBUST_fs_lights_hrs.dta", replace
}

****************************************************************** 
****************************************************************** 








